Collective excitations of trapped Bose condensates in the energy and time domains 
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A time-dependent method for calculating the collective excitation frequencies and densities of a 
trapped, inhomogeneous Bose-Einstein condensate with circulation is presented. The results are 
compared with time-independent solutions of the Bogoliubov-deGennes equations. The method is 
based on time-dependent linear-response theory combined with spectral analysis of moments of the 
excitation modes of interest. The technique is straightforward to apply, is extremely efficient in 
our implementation with parallel FFT methods, and produces highly accurate results. The method 
is suitable for general trap geometries, condensate flows and condensates permeated with vortex 
structures. 
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I. INTRODUCTION 

The equations denning the frequencies of single par- 
ticle excitations from a weakly-interacting dilute Bose 
as at zero temperature were introduced by Bogoliubov 
and applied to low-temperature Fermi systems by de 
Gennes pi . The underlying assumption is that tempera- 
tures are sufficiently low such that the condensate is not 
significantly depleted and that the interparticle interac- 
tions are weak and the gas dilute. Under such constraints 
the dynamics are dominated by single particle excitations 
meaning that semiclassical approximations can be used 
for the field excitations and the equations are equivalent 
to linear response theory g, ||, ||. Trapped hydrostatic 
ultracold gas clouds can be modeled to a high degree of 
precision under these conditions. Even within this frame- 
work the equations must be solved numerically because 
of the inhomogeneity and nonlinearity of the condensate. 
Determining the elementary excitations is much more 
complex in situations in which the condensate is flowing 
or permeated with defects such as vortex arrays. The 
extraordinary coldness of the atom cloud means the con- 
densate dynamics occur on millisecond time scales and 
thus are accessible to measurement as they evolve. Shape 
distortions, created by manipulating the condensate with 
external electromagnetic fields, can be used to measure 
mode frequencies || [?], |8| , wave-mixing || and damping 
rates In this way elementary excitation frequen- 

cies at low temperatures can be measured to an accuracy 
within a few percent, and the agreement with theory for 
the low-lying modes is astonishingly good [||. Precise 
measurements closer to the critical temperature are much 
more uncertain, and indeed theoretical modeling is also 
challenged by the influence of complex pair excitations 
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flllfl . The agreement between theory and experiment is 
much less satisfactory in this region. 



The experimental and theoretical study of condensates 
in the time domain mirrors experimental breakthroughs 
in ultrashort laser pulses used to probe ultrafast pro- 
cesses in chemistry and biology. Dynamic simulations 
and experiments applied to Bose condensates have pro- 
vided insight into the transition from superfluid flow to 
dissipation and drag [|l2| |l3) and in condensate forma- 
tion and destruction. Trapped condensates are partic- 
ularly suited to spectral methods, specifically the split- 
operator technique |L4j combined with the FFT method 
|15| . We report results in which efficient parallel comput- 
ing techniques have been implemented within the spec- 
tral method to determine the elementary excitation fre- 
quencies and densities. The method is ideally suited to 
the analysis of elementary excitations in complex flows 
and in the presence of defects. The method also allows 
the excitation mechanism of experiments to be modeled 
realistically and thus help devise schemes which produce 
optimal excitation of certain modes or superposition of 
modes leading to squeezed and entangled states jl7j , and 
observe collective excitations of quadrupoles giving evi- 
dence of superfluid motion Jl8], 



The paper is structured as follows: In section II the 
units and notation are described, and the field equations 
are introduced. The case of excitation in the presence 
of condensate circulation is discussed. The equations de- 
scribing excitations in a cylindrically symmetric trap and 
quadrupole excitations in an asymmetric trap are given. 
In section [II the time-dependent method is formulated. 
Finally, in sections IV and |v|, results are compared be- 
twen the two methods and some general conclusions are 
given on the relative merits of each method. 
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Interaction Strength 
C 


Chemical potential (p) 
K = K = 1 





1.000000 


2.000000 


10 


1.620434 


2.368791 


50 


3.020229 


3.489876 


100 


4.143006 


4.513850 


250 


6.415956 


6.685521 


500 


9.003106 


9.213175 


1000 


12.678320 


12.840724 



TABLE I: Values of the chemical potential (/i) of the k = 
and k = 1 states for the isotropic 2D trap as a function of 
interaction strength, C, in units of Tiloq 





qe = 


- —2 


qe ■- 


= 2 


c 


TDLR 


DVR 


TDLR 


DVR 





0.0000 


0.0000 


2.0000 


2.0000 


50 


0.8333 


0.8331 


1.8647 


1.8649 


100 


1.0299 


1.0300 


1.7752 


1.7748 


250 


1.1876 


1.1871 


1.6590 


1.6594 


500 


1.2565 


1.2560 


1.5902 


1.5906 


1000 


1.3027 


1.3022 


1.5388 


1.5393 



TABLE IV: Comparison of TDLR and DVR results for 
quadrupole excitation q r = 0,qg = ±2 for a condensate with 
circulation k = 1. 





K - 


= 




= 1 


c 


TDLR 


DVR 


TDLR 


DVR 





1.9995 


2.0000 


1.9996 


2.0000 


50 


1.9998 


2.0000 


1.9998 


2.0000 


100 


1.9936 


2.0000 


1.9979 


2.0000 


250 


1.9977 


2.0000 


1.9980 


2.0000 


500 


1.9978 


2.0000 


1.9978 


2.0000 


1000 


1.9978 


2.0000 


1.9977 


2.0000 



TABLE II: Comparison of TDLR and DVR results for the 
excitation frequencies uij, for the mode q r = l,qe = with 
k = 0, 1. 



low temperatures, treats the condensate semiclassically 
while quantizing the excitations so that: 

*(r) = VN <j)(r) +tp(r) (2) 

The condensate wavefunction satisfies the Gross- 
Pitaevskii equation: 

H Q (r)<P(r) + NU \^{r)\ 2 4>{r) = (3) 
Then expanding i\> in modes 

i>(r) = ]T [ujW&j + v*(r)a]\ (4) 



II. TIME-INDEPENDENT QUANTIZED FIELD 
EQUATIONS 



The dilute system of N bosonic atoms mass m is 
trapped by external fields 14 x t(r) and interacts weakly 
through the two-body V(r, r'). If the external field is the 
static trapping potential, the Hamiltonian for the system 
can be written 

H' = j dr *t(r)fT *(r) 

\ j J drdr'* t (r)* t (r')l / (r,r')*(r')*(r)(l) 



+ 



with Hq 



-V 2 + Vtrap(f ) — Mj with chemical potential 



/i. The interactions can be represented perturbatively 
by the pseudopotential V(r,r r ) = Uq5(i — r') where the 
interaction is proportional to the s-wave scattering length 
a, Uq = A-Kh 2 a/m. Bogoliubov's approximation, valid at 



c 


TDLR 


DVR 





2.0000 


2.0000 


50 


1.5569 


1.5569 


100 


1.5003 


1.5002 


250 


1.4561 


1.4561 


500 


1.4380 


1.4380 


1000 


1.4277 


1.4276 


oo 


1.4142 


1.4142 



TABLE III: Comparison of TDLR and DVR results for the 
quadrupole excitation q r — 0, qe = ±2 mode with condensate 
k = 0. 



The eigenvalue problem reduces to the Bogoliubov- 
deGennes equations in the form: 

£v,j(r) + NUo<j) 2 Vj(r) — hcjjUj(r) (5a) 
Cvj{r) + NU (j)* 2 u 3 {r) = -hu 3 v {r) (5b) 

where C = Hq(t) + 2NUo\<fr(r)\ 2 . Time-reversal sym- 
metry is reflected in the fact that every set of solutions 
{Ej,Uj,vj} has a corresponding set {—Ej,u*,v*}. For 
Ej > the functions Uj , Vj are orthogonal with normal- 
ization 



j dr(| Uj -(r)| 2 -K(r)| 2 ) 



(6) 



In our model the trap is represented by an asymmetric 
harmonic well: 

t r i \ 1 2/2 2, 22, 22\ 

V tlap (r) = ±muj (e x x + e y y +e z z) 

where e x ,y,z are restoring force strengths, and loq the nat- 
ural angular frequency. Numerical calculations are car- 
ried out in scaled dimcnsionless units. Length, time and 

energy are given in units : (Ti/2mu>o) 2 , uj and Hujo, 
respectively. 

When particle interactions are weak, corresponding to 
the ideal gas limit, the energies for an asymmetric trap 
are given by the oscillator formula 

E n = (n x + \)e x + (n y + \) e y + (n z + |) e z 

with n Xl n yi n z — 0,1,2,.... For a 2D symmetric os- 



cillator (e 2 



0, e x 



1) the cylindrical quantum 



numbers (q r ,qg) can be used: 

E n = n x + n y + 1 = 2q r 
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FIG. 1: Quadrupole moment evolution f xy (t) for low-lying 
excitations of a condensate with k = 0. 
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FIG. 2: Quadrupole moment evolution f xy (t) for low- lying 
excitations of a condensate with k = 1. The figure displays 
interference between the modes qg — ±2 



FIG. 3: Power spectral density of quadrupole moments (xy) 
versus angular frequency. The vertical scale is logarithmic 
and ranges from 0.1 to 1.0. Results for k = 0, dashed line, 
qg — ±2 and k = 1, solid lines, qg = ±2 for interaction 
strength C = 1000. 



where C — 8TraN a with N a the linear density of atoms 
in the z direction. The excitations can be written: 

U] (r,6) = u qr {r)e l{qe+K)0 (9) 
v 3 (r,9) = v qr (r)e i{qi> ~ K)e 

where j — {q r ,qg}- The eigensystem equations ( pajphj ) 
become 



d 2 l_g_ (gg_+ w) 2 
<9r 2 r dr r 2 



Xr 2 - n + 2C|0| S 



V( r ) 



+C(/> 2 S gr (r) =w 9r M ?r (r) (10) 



with g r = 0,1,2,... and qe = 0, ±1,±2, ... so that 
the excitation frequencies above the ground state (q r = 
0, qg = k) are: uij = 2q r + \qg + k\ — k. 



d 2 Id (qg - k) 



1^2 



r 2 -fi + 2C\(p\< 



dr 2 r dr r 2 4 

+C4> 2 u qr {r) = ~uj qr i qr {r) (11) 



B. Discete variable representation 



A. Excitations of vortex states 

A precise test is provided by low-dimensional models 
for which highly accurate results can be obtained. It also 
allows the study of excitations in the presence of vortices. 
For a single vortex line with circulation k along the axis 



of a cylindrically symmetric trap e x = e y = 1, e 2 



0: 



(7) 



Equations (|8|),(|l0|) and ( |ll| ) are solved using the dis- 
crete variable representation (DVR) 20 1 . The condensate 
and excitation functions are represented by a Lagrange 



mesh (j>, u qr , v qr . 



N 



i=i 



N 



,(r) = ^u^r^f^r) 



(12) 



(13) 



d 2 1 d 
dr 2 r dr 



(8) 



C|0| 2 0(r) = Mi?) 



N 



(14) 
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where fi(r) are Lagrange interpolating functions con- 
structed from a set of orthonormal functions 



Xk(r) = h k 2 w{r)ip k {r) 



such that 



N-l 



Mr) = xl(n)xk(r) 



(15) 



(16) 



We employ generalized Laguerre polynomials J2jj ] 
p k (r) = L%(r) = ((a + k)\/k\a\)iFi(-k;a + l;r) asso- 
ciated with the weight function w(r) = r a e~ r and nor- 
malization factor hu = [k + oi)\/k\. The representation of 
the centrifugal energy terms of equations (fl0|) , (|ll|) can be 
marginally improved by using different values of a [ p0[ . 
Different grids for each angular momentum state requires 
the use of interpolation to connect different states for 
k =/= 0. However using a common grid is more practical. 
The DVR results presented correspond to k = N = 50 
and the common value a — 1 so that: 



Mr) 



1 XN(r) 



X N ( r i) 



(17) 



where the mesh points r, are the N zeros of L%(r). Ap- 
plied to the coupled equations (|§|) , (|To|) , (pT]) this method 
reduces the problem to that of solving a system of nonlin- 
ear equations for the condensate and a linear eigenvalue 
problem for the excitations. The set of (N + 1) nonlin- 
ear equations is efficiently solved using Newton's method 
where \x is an unknown and normalization of cj> is im- 
posed. As the functions 0, u, v are to be represented on 
the same grid points, an ade quat e coverage of points 
is required. A scaling factor h |20) is used to contract 
the natural mesh so that it extends to ~ 1.5rxF where 
t*tf ~ 2(C/27r) 1 / 4 is the Thomas-Fermi radius. Increas- 
ing the number of mesh points beyond N — 50 does not 
improve the accuracy beyond the figures quoted. While 
the matrices are dense they give results far more accurate 
than finite difference methods using similar sized matri- 
ces. We used a standard library routine for solution of 
the eigenvalue problem |24|. 



III. TIME-DEPENDENT LINEAR RESPONSE 
EQUATIONS 

The mean-field approximation to the field operator is 
ip(r,t) « (>J>(r,i)) w (^(r,t)). The angular brackets 
denote averaging with respect to highly-occupied (N ^> 
1) condensate number states. The evolution equation 
for the mean field is the time-dependent Gross-Pitaevskii 
equation |J: 



The external potential includes the trapping potential 
and a small time-dependent perturbation: V cx t{r,t) — 
Vt iap (r) + Vpert (r,t). The system responds to this per- 
turbation by populating modes selected by the symmetry 
of the perturbation: 

i/>(r,t) « e- ifit (p(r) (19) 

+ K'W^(r)e-^*- vt +b J (t)v J (r)e lut -^ t ] 
i 

The oscillations of the system provides data for the fre- 
quencies: ujj. This replaces the nonlinear coupled equa- 
tion eigenvalue problem by an initial value problem fol- 
lowed by spectral analysis. Alternatively, we can intro- 
duce a perturbed initial state ^(r, 0) seeded with the 
spatial form to overlap with the symmetry of the de- 
sired excited state. This is equivalent to a sudden per- 
turbation of the system. To study low-energy excitations 
the ground state, <f>(r), is found by propagation of equa- 
tion ( |l8|) in imaginary time fl3|| using a trial wavefunc- 
tion. Implementing the split-operator FFT method on 
a parallel computer using between 1 and 16 processors. 
The method is extremely well-suited for nonlinear evolu- 
tion equations |TJ| . In imaginary time the highly-excited 
modes are exponentially suppressed, and while it is not 
necessary to eliminate these modes before introducing 
the perturbation, a pure ground state is more efficient 
for spectral analysis. 

The dipole modes were excited by a sudden but small 
displacement of the trap center. The subsequent center- 
of-mass oscillations are known exactly (tu = e Xl e y ,e z ) 
and should provide precise benchmarks for the numer- 
ical method employed. Breathing (monopole) modes 
were excited selectively, and in line with experimental 
methods ||, by adjusting the force constants to com- 
press or expand the condensate e X)H — ► 1 ± e. The 
low-lying quadrupole modes are accessible by seeding 
the initial state with a term of the appropriate sym- 
metry. For example, adjusting the true ground state to 
^(t-jO) = <j>(r)(l + exy), and following the evolution of 
fy(r,t) under the Gross-Pitaevskii equation, the beats 
of the xy quadrupole spectrum are evident. While sin- 
gle modes may be selected by adiabatic perturbations 
of selected frequencies, sudden and strong perturbations 
can be used to gain data on higher excited modes. A 
significant advantage of the method is that a variety of 
perturbations of different symmetries can be used simul- 
taneously to obtain the excitations by a single evolution 
of the nonlinear equations. 

The wave packet method combined with spectral 
method was pioneered by Feit and Fleck and has 
been used widely |l6|, ^2], 23 . The response of the sys- 
tem is measured by the moment: 



_ d , , . 



h 

_ v 2 + y cxt (r,t)+A[/ |^(r,i)| i 
2m 



ip(r,t) 
(18) 



f A (t) - J drV*(r,t)A(r)y(r,t) 



(20) 



The excitation frequencies emerge from the Fourier trans- 
form of /a(0- The finite sampling time means that the 
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power spectral density signal can be improved by win- 
dowing. The windowing function we use is a /3-valued 
Kaiser function chosen to minimize the background. The 
local maximum of the power spectrum is then used to 
identify the frequencies. The density profiles of the quasi- 
particles can be extracted by further spectral analysis. 
Once the appropriate uij is known the Fourier transform: 

gf{r) = T" 1 / T e t ^ ± ^)**(r,i) dt isolates the corre- 
sponding amplitude. For both the imaginary time and 
real time propagation we used the same split-operator 
scheme. In the calculations we used a grid of 64 points 
in each dimension. To create the ground state we used 
a variable time step along with extrapolation to increase 
the speed of calculation and control the accuracy of the 
ground state. To obtain an accurate spectral resolution 
a total propagation time of ~30 times the trap period 
is necessary. In addition, the signal was processed with 
a windowing mask. Alternatively if the perturbation is 
small and introduced slowly, on the timescale of several 
trap periods, only a very small number of modes will be 
excited. The correlation /a(£) can be fitted to a Fourier 
expansion in just a few frequencies. This has been done 
for all the frequencies presented in this paper and gives 
equivalent accuracy while requiring only 3 trap periods. 



A. Classical limit 

When interparticle forces dominate the quantum pres- 
sure (C — » oo) and the hydrodynamic limit is reached. 
The equilibrium condensate density, p = \ip\ 2 , for n = 
corresponds to the Thomas-Fcrmi distribution: 



Po (r) = (p-V(r))/C 

Acoustic modes of excitation [^5] p(r, t) = 
p'(r) e~ tuJt are determined by the equation 

u?p' = -2CV • (poVp') 



(21) 
Po(r) + 

(22) 



For the 2D cylindrical trap (e x = e y = l,e z = 0) this 
gives (C — > oo): 



wj? = 2g r (? r + l) + g e (2? r + l). 



(23) 



where the radial and angular quantum numbers (q r ,qe) 
are as defined above. 

For a 3D asymmetric trap (e x ^= e y ^= e^) the lowest 
quadrupole modes of scissors motion H H that is 
of the form p' cx xy, yz, xz are given by the solutions: 



^xy e x e y ^yz 



el + el, andwL 



IV. RESULTS 

In this article we compare the time-dependent (TDLR) 
and time-independent methods (DVR) for excitation in 
2D traps. For the DVR calculations we used N = 30 
and N = 50 radial mesh points with a scaling factor 



c 






LOyz 





3.0005 


4.0003 


5.0002 


10 


2.8059 


3.8041 


4.7938 


50 


2.5439 


3.5536 


4.4612 


100 


2.4432 


3.4228 


4.2521 


250 


2.3521 


3.3134 


4.0254 


500 


2.3101 


3.2609 


3.8951 


1000 


2.2821 


3.2257 


3.7993 


oo 


2.2361 


3.1623 


3.6056 



TABLE V: Lowest-energy quadrupole modes for an asym- 
metric trap with e x = 1, e y =2 and e z = 3 as a function of 
interaction strength C. Results for C = oo correspond to the 
hydrodynamic limit 



h = 0.08. The results presented, using N = 50, are con- 
verged to at least 6 decimal places. For reference, the 
values of the chemical potential are presented in Table 
[j]. The TDLR calculations were performed on grids of 
64 x 64 x 64 points. Convergence and accuracy can be 
calibrated by comparison with exact results. For exam- 
ple, the 2D dipole (q r = 0,qg = 1) mode is a feature 
of the trap and independent of the particle interactions 
(C) and circulation k, so that uuj = 1. The DVR results 
for this mode are accurate to at least six decimal places, 
while the TDLR results are within 0.3% of the exact val- 
ues over the range < C < 1000. Table shows the 
results for the lowest breathing mode (q r = 1, qg = 0) for 
k = and k = 1. Again, in 2D this mode is not sensi- 
tive to particle interactions. TDLR results are stable and 
consistent for this mode as C increases, slightly under- 
estimating the frequency in the fourth significant figure. 
The loss of precision in the TDLR results for high C is 
of the same order as the errors in the chemical potential 
results. 

The quadrupole modes are more sensitive to changes 
in C and also to the condensate circulation. In figure 
[l] the correlation f X y{t) is shown for C — 1000. The 
smooth sinusoidal curve translates to the Fourier trans- 
form spectral density shown as the dashed line in fig- 
ure [|. The spectral density shows a noise-free profile 
with a well-defined peak at the value u> — 1.4277. As 
C — > +oo, lo — > \/2. In Table [II results are presented 
for the mode: k = 0, q r = 0, qg = ±2 . The agreement be- 
tween the time-dependent and time independent results 
is extremely good over the entire range of C. Given that 
the dipole mode showed the TDLR accuracy is within 
four figures, the accord to five figures for C = 250 and 
C = 500 is fortuitous. 



When condensate circulation is present the quadrupole 
circulation can flow with or against the background and 
the qg = ±2 modes split. This is evident in the beats aris- 
ing in the f X y(t) correlation figure |^ and the correspond- 
ing spectrum (figure ||[). The agreement with DVR re- 
sults, indicated in Table IV is excellent. Since the TDLR 



method is based upon a cartesian grid, it is readily used 
for arbitrary geometries of the trap. For example, the 
quadrupole modes xy,yz, zx of an asymmetric 3D ellip- 
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soidal trap can found in the same way. In this case we 
used a mixed perturbation of all three symmetries and 
then frequency analyzed the correlations of each moment. 
The results are presented in Table 0. The data tends to 
the correct high C limit uniformly and accurately. 

V. CONCLUSION 

We have compared methods for calculation of the ex- 
citation frequencies of a condensate at zero temperature. 
The method based on linear response theory is extremely 
efficient and successful at producing accurate results and 



can be tailored to the experimental methods used to cre- 
ate the excitations of interest. The method is straight- 
forward to implement, yields reliable results, and can be 
computed cheaply and efficiently. The versatility of the 
method extends to complex trap geometries and conden- 
sate topologies. Thus it is suitable for the study of ex- 
citations and sound propagation in traps where the con- 
densate contains soliton and vortex structures |p7f . 
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tational Physics Laboratory in the framework of the Eu- 
ropean Community - Access to Research Infrastructure 
action of the Improving Human Potential Programme. 
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